'there should be data already set up as series for b, p, q, pop, mu, gini, theil for each country
smpl 1 89
genr lamda=pop/@sum(pop)
'scalar test=@sum(lamda)
genr muuu=lamda*mu
scalar muu=@sum(muuu)
genr s=lamda*mu/muu
genr lsg=lamda*s*gini
genr l22=lamda*log(b)
scalar l2=@sum(l22)
genr l33=lamda*(@digamma(q)-@digamma(p))
scalar l3=@sum(l33)
genr lw11=lamda*log(mu/b)
scalar lw1=@sum(lw11)
genr lb22=lamda*log(mu)
scalar lb2=@sum(lb22)
'calculate Theil's L index
scalar l=log(muu)-l2+l3
scalar lw=lw1+l3
scalar lb=log(muu)-lb2
'the following modified from Bill 's program
genr p1=p+1
genr q1=q-1
scalar gb=0
scalar gw=0
scalar slum=0
smpl 1 50000
for !i=1 to 89
scalar pp=p(!i)
scalar qq=q(!i)
genr x=@rbeta(pp,qq)
scalar bb=b(!i)
genr y=bb*x/(1-x)
for !j=1 to 89
scalar ppp=p(!j)
scalar qqq=q(!j)
scalar bbb=b(!j)
genr xx=y/(y+bbb)
genr fy=y*@cbeta(xx,ppp,qqq)
scalar m=@mean(fy)
scalar slum=slum+lamda(!j)*lamda(!i)*m
scalar gb=gb+lamda(!j)*lamda(!i)*@abs(mu(!j)-mu(!i))
next
next
'calculate regional gini total, gini within, between and the left over
scalar g=-1+(2/muu)*slum
scalar gw=@sum(lsg)
scalar gb=gb*(1/(2*muu))
scalar gr=g-gw-gb

